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Models of two-dimensional (2D) traps, with the double-well structure in the third direction, for 
Bose-Einstein condensate (BEC) are introduced, with attractive or repulsive interactions between 
atoms. The models are based on systems of linearly coupled 2D Gross-Pitaevskii equations (GPEs), 
where the coupling accounts for tunneling between the wells. Each well carries an optical lattice 
(OL) (stable 2D solitons cannot exist without OLs). The linear coupling splits each finite bandgap 

| in the spectrum of the single-component model into two subgaps. The main subject of the work is 

£^ , spontaneous symmetry breaking (SSB) in two-component 2D solitons and localized vortices (SSB 

was not considered before in 2D settings). Using variational approximation (VA) and numerical 

' methods, we demonstrate that, in the system with attraction or repulsion, SSB occurs in families 

of symmetric or antisymmetric solitons (or vortices), respectively. The corresponding bifurcation 
destabilizes the original solution branch and gives rise to a stable branch of asymmetric solitons or 

' vortices. The VA provides for an accurate description of the emerging branch of asymmetric solitons. 

| In the model with attraction, all stable branches eventually terminate due to the onset of collapse. 

. Stable asymmetric solitons in higher finite bandgaps, and vortices with a multiple topological charge 

1 are found too. The models also give rise to first examples of embedded solitons and embedded vortices 

(the states located inside Bloch bands) in two dimensions. In the linearly-coupled system with 
opposite signs of the nonlinearity in the two cores, two distinct types of stable solitons and vortices 
are found, dominated by either the self-attractive component or the self-repulsive one. In the system 
Mh ' with a mismatch between the two OLs, a pseudo-bifurcation is found: when the mismatch attains its 

£H , largest value (tt), the bifurcation does not happen, as branches of different solutions asymptotically 

approach each other but fail to merge. 

PACS numbers: 03.75.Lm, 05.45.Yv, 42.65.Tg, 47.20.Ky 
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I. INTRODUCTION 

00 , 
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■ A great variety of dynamical patterns in Bose-Einstein condensates (BECs) can be supported by optical 
lattices (OLs), i.e., spatially periodic potentials induced by laser beams illuminating the BEC [1]. A well- 

qq . known species of such patterns are gap solitons in condensates composed of repulsively interacting atoms. 

■ Gap solitons in BEC were predicted theoretically 0-0, and then created in a "cigar-shaped" (i.e., nearly one- 
! dimensional, ID) crossed-beam optical trap, to which a longitudinal OL was added. In self-attractive media, 

the OL potential makes it possible to trap a usual soliton at a required position, and generates stable multi- 
soliton complexes In the case of a very deep OL, the corresponding Gross-Pitaevskii equation (GPE) 0] 
reduces to a discrete nonlinear Schrodinger (NLS) equation [8], which supports well-known solitons, including 
/\ ' staggered ones 0|. 

Another topic which has drawn a great deal of interest, in terms of BEC [10] and nonlinear optics 
alike, is spontaneous symmetry breaking (SSB) in matter waves or optical beams trapped in settings based 
on double- well potentials (using the terminology widely adopted in optics the medium corresponding 
to each potential well is sometimes called a "core" below). If two parallel "cigar-shaped" traps are strong 
enough, the system may be described by a pair of one-dimensional GPEs with linear coupling between 
them [12], hence the destabilization of symmetric solitons and spontaneous transition to asymmetric ones 
in the double cigar-shaped trap, filled with self-attractive condensate, is described in the same way as the 
formation of asymmetric solitons in the well-studied model of dual-core nonlinear optical fibers [l6[ (hi 
nonlinear optics, the SSB of ID solitons was also studied in dual-core models with non-Kerr nonlinearities, 
such as quadratic [17] and cubic-quintic fl8jl). On the other hand, if the dual-core BEC waveguide is realized 
as a set of two parallel finite- width potential troughs, separated by a finite buffer layer, which are embedded 
in a 2D (horizontal) trap and tightly confined in the transverse (vertical) direction, the SSB in the respective 
double-core solitons can be examined in the framework of the full 2D GPE (i.e., considering the tunneling 
across the buffer layer explicitly, rather than approximating it by the linear coupling between the two ID 
equations) [l9[. Both models predict a subcritical symmetry-breaking bifurcation for solitons in the self- 



2 



attractive condensate [which means that branches of asymmetric solutions emerge as unstable ones at the 
SSB point, and originally go backward (for which reason the bifurcation is also called a backward one), but 
then turn forward, stabilizing themselves at the turning point]. 

The objective of the present work is to find symmetric, antisymmetric and asymmetric families of 2D BEC 
solitons, and similar families of localized vortices (alias vortex solitons), in the system built as a stacked 
pair of "pancake-shaped" traps, each carrying a two-dimensional OL. The traps are linearly coupled by 
tunneling across the buffer layer separating them. This model is a natural extension of a double-trap system 
recently studied in the ID setting [12[ (in optics, spontaneous symmetry breaking in ID multi-component gap 
solitons was studied in models of dual-core [l3[ and triple-core [14[ fiber Bragg gratings). The generalization 
to the higher dimension is significant, because 2D solitons are drastically different from their ID counterparts 
[l5[, especially in the case of the self-attractive nonlinearity (the existence of the respective 2D solitons is 
fundamentally restricted by the possibility of collapse). Besides that, the new species of vortex solitons is 
possible in 2D geometry |15j . To the best of our knowledge, SSB in 2D solitons or vortices has never been 
considered before. 

We will examine different physically relevant versions of the model, with attractive or repulsive nonlinearity 
in each trap (following Ref. [12[, they are referred to as AA and RR systems, i.e., "attractive-attractive" 
and "repulsive-repulsive"). In the AA system, the solitons originate in the semi-infinite gap of the linear 
spectrum, while in the RR system one may expect to find solitons in finite bandgaps induced by the OL 
(we will consider the first two gaps, each of them being split into two subgaps under the action of the linear 
coupling between the traps). In some cases (similar to what was found in the ID model [12[), families of 
solitons and vortices extend across Bloch bands separating the (sub)gaps, thus becoming embedded solitons 
po| . As a matter of fact, the present paper reports the first examples of 2D embedded solitons (and vortices). 

A mixed (AR) system, with self-attraction in one trap and self-repulsion in the parallel one, will be 
considered too. As is known, the sign of the interaction between atoms can be reversed by means of the 
Feshbach resonance imposed by external magnetic field [21]; accordingly, the signs may be made opposite in 
the two traps by applying the magnetic field which is inhomogeneous across the "pancake" stack. Solitons 
in ID settings of the latter type, but without OLs, were elaborated in Refs. (22[; a 2D generalization was 
considered too (with no OLs either) but it does not give rise to stable solitons. 

In the normalized form, 2D models considered in this work are based on a system of linearly coupled GPEs 
for the mean-field wave functions in the two flat traps, ip(x,y,t) and y, £), 

i^ + V 2 ^ + ^[cos(2x) + cos(2?/)]^ + Ai|^| 2 V ; + ^ = 0, 

(la) 

i(j>t + V 2 </> + £ [cos(2:r + Ai) + cos(2?/ + A 2 )]</> + A 2 |(/>| 2 (/> + = 0, 

where x and y are planar coordinates, V 2 = d 2 /dx 2 + d 2 /dy 2 , Ai 5 2 = +1 and —1 correspond to the 
attractive and repulsive nonlinearity in the cores, real linear-coupling coefficient n is proportional to the rate 
of tunneling across the potential barrier separating the parallel traps ("cores"), and e is the strength of the 
OLs, whose period is scaled to be tt. We assume equal depths of the OLs in both traps, while asymmetry 
between the lattices may be introduced by mismatches (phase shifts) in the x- and ^/-directions, Ai and A 2 . 
We will chiefly deal with the symmetric system, i.e., Ai = A 2 = 0; nevertheless, effects produced by the 
mismatches will be addressed too (in fact, we will consider the case of the largest possible mismatch in the 
diagonal or horizontal direction, i.e., Ai = A 2 = 7r, or Ai = 7r, A 2 = 0, respectively). In optical models, 
various effects of the mismatch between linearly-coupled Bragg gratings on gap solitons supported by the 
gratings were studied in Refs. [24j ]. 

The 2D form of Eqs. ([Taj) implies that the "pancake" -shaped traps are tightly confined in the transverse 
(Z) direction, which allows one to reduce the underlying 3D GPE to two dimensions, as performed, in the 
general form, in Refs. (25[. Together with normalizations of the OL period and nonlinearity constants 
adopted in Eqs. ([Taj) , this implies that variables t and (x, y) in Eqs. ([Taj) are related to physical time T and 
coordinates (X,Y) by 

^2^ T ' (2) 
where m is the atomic mass, and d the OL period in physical units, while the scaled wave functions are 
related to their counterparts, ^ and <£, in the 3D space as follows: 



Y,Z,T)\_1_ f iP{x, t) l (_L wlT - ^zA 

*(X,Y,Z,T) / " 2V \a s \d? \ 4>{x,t) / 6XP 1, 2^ 2h Z ) 



(3) 
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Here uo±_ is the transverse trapping frequency, and a s the s-wave scattering length of atomic collisions (a s < 
for attractive interactions). Due to the normalizations, the lattice strength is represented by e = 2£/£ rec , 
where £ rec = (ttH) 2 / (md 2 ) is the lattice recoil energy, and £ the depth of the periodic potential in physical 
units. 

The A A and RR systems may be experimentally realized in 7 Li [26] or 85 Rb J2?J, and 87 Rb [1] condensates, 
respectively. In this work, the analysis is presented in the range of ft ~ 1 in Eqs. ([Taj) . With regard to 
Eqs. ([2]), the corresponding normalized tunnel-coupling time, t couv> \ = 7r/(2ft), translates, for d ~ 1 /im, into 
physical time T coup i ~ 10 fis and 100 /is, for lithium and rubidium, respectively. A prediction most essential 
to the experiment is the number of atoms expected in stable solitons. Assuming an experimentally relevant 
value of the trapping frequency - for instance, uo±_ ~ 2tt x 100 Hz for the rubidium condensate (then, the 
corresponding transverse-confinement size is estimated to be ~ 1 /im, i.e., comparable to the OL period) - 
and translating the results reported below into physical units, we conclude that the symmetry-breaking 2D 
solitons are built of ~ 10 3 atoms. In vortex solitons, the number of atoms is larger, roughly, by a factor of 
10. 

The paper is organized as follows. In Section II, we focus on relevant mathematical formulations, which 
include the linear spectrum of the coupled system, stationary equations for solitons, and the setting for the 
variational approximation (VA) describing symmetric and asymmetric solitons in the A A and RR systems. 
Systematic numerical results for solitons in symmetric AA and RR systems (and comparison with predictions 
of the VA) are reported in Section III, while Section IV reports numerical results for vortex solitons, in the 
same systems. Findings for asymmetric systems (of the AR type, as well as ones with mismatched OLs) are 
collected in Section V. Section VI concludes the paper. 

II. FORMULATIONS: THE LINEAR SPECTRUM, STATIONARY EQUATIONS, AND 

VARIATIONAL APPROXIMATION 

A. Linear spectra 

Before looking for soliton solutions, it is necessary to identify the system's spectrum within the framework 
of the linearized equations. First, we consider the symmetric system, with Ai = A2 = 0, and look for 
solutions with chemical potential ji as 

{^(z, y, t), 0(£, V, t)} = [a(x, y) ± /3(x, y)] e~^. (4) 

Substituting this in the linearized version of Eqs. ([Taj) . we arrive at a 2D eigenvalue problem based on the 
decoupled equations, 

V 2 a(x, y) + e [cos(2x) + cos(2?/)] a(x, y) = — (fjL + K,)a(x,y), (5a) 
X7 2 [3(x,y) +s[cos(2x) +cos(2y)]P(x,y) = -(jjl - k)/3(x, y). (5b) 

Each one of these equations is tantamount to the eigenvalue problem for the single GPE with a 2D lattice and 
effective chemical potential \i' = \i±k. The spectrum generated by the linearization of the single-component 
GPE is well known, see Refs. Q and references therein. Figure [1] shows a typical example of the spectrum of 
the 2D coupled system with the zero mismatch. We observe that the linear coupling in Eqs. ([Taj) splits the 
finite gaps of the single-component model BEC into pairs of subgaps, which is similar to what was observed, 
under the action of the linear coupling, in the ID counterpart of the present system [12|. 

To show in detail how the gaps split into subgaps, in Fig. [2] we display the dependence of the chemical 
potential on vectorial quasi- momentum (k x ,k y ) in the diagonal direction (k x = k y = fc), in the uncoupled 
system and in its coupled counterpart with ft — 1, both pertaining to e = 8. The new, nearly flat, branch of 
the fi(k) dependence, which appears in the latter case around \i = —7.5, explains the splitting of former gap 
1 into subgaps la and lb by a new very narrow Bloch band, cf. Fig. [TJ At other values of £, the splitting 
mechanism is quite similar. 

Figure [2] and all others in this paper display numerical results for OL strength e = 8 (a moderately deep 
lattice), which corresponds to the most typical situation. As for coupling constant ft, the results are displayed 
for ft = 1 (which also helps to presents generic findings), unless a different value of ft is indicated. 

As shown in Fig. O we also examined the transformation of the linear spectrum under the action of 
the diagonal mismatch between the OLs in the two traps, assuming Ai = A2 = A. In this case, the 
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FIG. 1: (Color online) The spectrum of the linearized 2D system with aligned optical lattices 
(Ai = A 2 = 0). Throughout the paper, all numerically obtained results are given for e = 8, which makes it 

possible to display the findings in a generic form. Here and in figures following below, Bloch bands are 
shaded, (a) The uncoupled system, n = 0; (b) the coupled one, with n = 1 (unless specified otherwise, all 

figures in the paper pertain to k = 1). 
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FIG. 2: (Color online) The dependence of chemical potential fi on quasi- momentum k (in the diagonal 
direction, see text) in the uncoupled system, n = (a), and in the coupled one, n = 1 (b). In both cases, 
Ai = A2 = and e = 8. The new nearly flat branch of the dependence, observed around \i = —7.5, 
explains the splitting of gap 1 into subgaps la and lb in Fig. [lj 



linearized equations cannot be decoupled; nevertheless, the computation is straightforward, as we can use the 
separability of the linear operator to construct eigenfunctions of the 2D model as products of eigenfunctions 
of the corresponding ID models [12[. Figure [3] demonstrates that, with the increase of A, the subgaps 
generated by the coupling-induced splitting tend to shrink, which resembles the trend observed in the ID 
model 0. 



FIG. 3: (Color online) The change of the linear spectrum with the increase of diagonal mismatch A 
between the lattices in the coupled traps, for n = 2. 



B. Stationary equations 

Stationary solutions to the full nonlinear system of Eqs. ([Taj) are looked for as {ip, <f\ = 
{u(x, y), v(x, y)} e~ Zflt , with real functions u and v to be found from equations 

\±u + V 2 u + e [cos(2x) + cos(2y)] u + Ai^ 3 + kv = 0, 

(6a) 

fiv + V 2 v + e [cos(2x + Ai) + cos(2?/ + A 2 )] v + X 2 v 3 + ku = 0. 

In the symmetric models, i.e., AA and RR ones (Ai = A2 = ±1) with zero mismatch (Ai = A2 = 0), 
symmetric (u = v) and antisymmetric (u = — v) solutions can be expressed in terms of a stationary solution 
of the single-component GPE with chemical potential /i, to be denoted as uo(x,y; ja): 

u(x, y; /i) = ±v(x, y\ fi) = u (x, y\v± «). (7) 

Similar to the situation in the ID model [l2[, we conclude from Eq. ([7j) that, when the gap splitting occurs, 
the symmetric and antisymmetric solitons will be moved to the lower and upper subgaps, respectively. In 
some cases, they may end up in Bloch bands separating the subgaps, thus becoming embedded solitons, 
examples of which are presented below. 

To find general asymmetric soliton solutions and SSB bifurcations linking them to their symmetric and 
antisymmetric counterparts, we solved the full system of stationary equations ([6a]) numerically. The stability 
of all solitons was examined by direct simulations of Eqs. ([Taj) . Soliton families are characterized (see below) 
by norms of the two components and the total norm, 

{N U ,N V } = J J {u 2 (x,y),v 2 (x,y)}dxdy, N = N U + N V , (8) 

For asymmetric solutions, we define the asymmetry ratio [l2[ as 

Q=\N U -N V \/(N U +N V ). (9) 
The total norm, together with the energy, are dynamical invariants of Eqs. ([Taj) . 

C. Variational approximation for the symmetric system 

The symmetric version of Eqs. ([6a]) (Ai = A2 = A, Ai = A2 = 0) can be derived from the following 
Lagrangian: 
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L = I I [fi (\u\ 2 + \v\ 2 ) - (\Vu\ 2 + \Vv\ 2 ) + e [cos(2z) + cospy)] (\u\ 2 + \v\ 2 ) 

+ (A/2) (M 4 + \v\ 4 ) + K (u*v + uv*)] . (10) 

A simple real isotropic Gaussian ansatz [28] , with single width W but different amplitudes A and B pertaining 
to the two components, may be adopted for the solitons: 

/ 2 2 \ 

{u,v) = {A,B)e^(- X -^-y (11) 

The total norm of this expression is [see Eq. ([8])] 

N = ttW 2 (A 2 + B 2 ) . (12) 
The substitution of ansatz (jTTJ) in Lagrangian (fTOj) yields the corresponding effective Lagrangian, 

LeffA = (A 2 + B 2 ) W 2 - {A 2 + B 2 ) + 2s (A 2 + B 2 ) W 2 e~ w " (13) 
+ (A/4) (A i + 5 4 ) W 2 + 2kABW 2 . (14) 

Variational equations following from this Lagrangian can be conveniently written as 

dL eS /d (A 2 + B 2 ) = dL ee /d (AB) = dL ee /d (W 2 ) = 0. (15) 

For symmetric and antisymmetric solitons, defined by A = sB with s = ±1, Eqs. ([T5]) amount to a set of 
two equations, 

TV 2ttW 2 A 2 =87r\(l-2sW /L e- w2 ^ , (16) 

fi + sk = -W~ 2 - 2ee~ w2 (l - 2W 2 ) . (17) 

In the attractive model (A = +1), Eqs. (fT6j) and (jTYJ) are tantamount to those analyzed in the framework 
of the VA in Refs. [28[ and (29[. It was demonstrated that the norm of the 2D soliton takes values in the 
following intervals: 

< TV < = 8tt, if e > e cr = e 2 /8 w 0.92; 

(18) 

N^ r) ee 8tt (1 - e/e cr ) < N < Sir, if e < e cr . 

Note that N^lx^ = 87r corresponds, in terms of the VA [30l |. to the norm of the Townes soliton, i.e., the 
localized solution of the 2D NLS equation in free space, which is unstable against collapse [3l| (while the 
OL stabilizes 2D solitons [H, H^). In fact, the vanishing of at e > e cr in Eqs. (fT8|) is an artifact 

of the VA [28], explained by the fact that Gaussian ansatz ([TT]) is not appropriate for multi-peaked solitons 
supported by the strong OL. 

In the repulsive model (A = —1), variational equations ([T6]) and (fT7|) which, essentially, pertain to the 
single-component setting, predict solitons only if the OL is strong enough [29[. Indeed, a straightforward 
consideration of Eq. ([T6]) demonstrates that, in this case, solutions exist only for e > £ cr , the norm of the 
solution family being limited to interval 

< JV < JV^> = 87T - 1) (19) 

[e cr is the same as in Eq. ([T8]) ]. The comparison of Eq. (fT9|) with known numerical results for 2D gap solitons 

0, [34| demonstrates that the nonexistence of solitons in the repulsive model at e < e& is another artifact 
of the VA, signalling that gap solitons cannot be approximated by simple Gaussian ansatz ([TT]) for e < £ cr . 
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For asymmetric solitons, with A 2 ^ B 2 , Eqs. (fT5|) can be cast in the following form: 



2 A( 



W 



-2 



2sW 2 e 




+ \l (W- 2 - 2eW 2 e~ w2 ) 2 + 2k 2 , 



(20) 



AB = 2Xk, 



(21) 



li = -2ee~ w2 (1 - W 2 ) - \\J (W- 2 - 2eW 2 e~ w2 ) 2 + 2k 2 . 



(22) 



With regard to Eq. (|H]) 



, norm (JT2J) of the asymmetric solitons becomes 




(23) 



Variational solutions for the asymmetric solitons are meaningful if they satisfy an obvious condition, 
A 2 + B 2 > 2 \AB\. In fact, the asymmetric solitons bifurcate from symmetric or antisymmetric ones [if, 
respectively, A = +1 or —1, as follows from Eq. ([2T]) ] precisely at point A 2 + B 2 = 2 \ AB\. Soliton families 
predicted by the VA, i.e., obtained by numerical solution of Eqs. ([2Q|) - ([23|h are represented by the corre- 
sponding dependences N = N(/i) in Figs. [Hand [71 along with their counterparts found from the numerical 
solution of full stationary equations (jSaj). 



In Fig. [H we display a generic example of families of antisymmetric, symmetric and asymmetric stationary 
solitons in the AA model with zero mismatch between the lattices (A = 0). The families were found from 
systematic numerical solutions of Eqs. (|6a|). Figure [4] also includes the solution families as predicted, for the 
same case, by the VA, which demonstrates good agreement between the variational and numerical results, 
for all the three types of solitons (at small values of the norm, the variational branches are indistinguishable 
from their numerical counterparts). A typical example of comparison of the profiles of asymmetric solitons 
produced by the numerical solution and VA is presented in Fig. [5l 

We observe in Fig. [41(a) that the symmetric-soliton branch undergoes a supercritical bifurcation at some 
critical value of the norm, giving rise to the branch of asymmetric solitons, which is the manifestation of 
the SSB in this setting (note that the bifurcation point is very accurately predicted by the VA). Symmetric 
solitons are stable before the bifurcation, and unstable afterwards, while asymmetric solitons emerge as 
stable solutions [recall the stability of the solutions was identified by direct simulations of Eqs. ([Taj) ]. 
On the other hand, the family of antisymmetric solitons never bifurcates, i.e., it never gets destabilized 
through SSB. As a consequence of that, the system features bist ability, when the antisymmetric solutions 
are stable simultaneously with either symmetric or asymmetric ones (below or above the bifurcation point, 
respectively). It is worthy to note that all branches of the solutions satisfy the Vakhitov-Kolokolov (VK) 
criterion, dN/djj, < 0, which is a necessary stability condition [3l|, [32[. It is also noteworthy that the stable 
branch of antisymmetric solitons displayed in Fig. Ufa) continues across the (narrow) Bloch band separating 
the semi-infinite gap and subgap la. Inside the narrow band, this family represents 2D embedded solitons 
(which are stable in direct simulations). As far as we know, this is the first example of 2D embedded solitons 
reported in any model. 

The above results resemble findings recently reported in the ID system [l2j . However, in the 2D case we 
observe an additional destabilization mechanism in the AA system, through collapse of the soliton of any 
type, when its norm becomes too large. In particular, the antisymmetric branch, which would be totally 
stable in the ID case, loses its stability in the 2D system when the soliton's norm exceeds the corresponding 
collapse threshold. Direct simulations demonstrate that, in this situation, the unstable antisymmetric soliton 
is first transformed into a strongly asymmetric structure; eventually, the high-amplitude component collapses, 
forming a singularity, while its counterpart with the lower amplitude decays into quasi-linear waves, as shown 
in Fig. [H 

The evolution of symmetric solitons destabilized by the SSB bifurcation is not affected by the collapse. 
Similar to what was reported in the study of the ID model [12[, the unstable symmetric solitons clearly tend 
to rearrange themselves into stable asymmetric counterparts (examples are not shown here, as they are not 



III. NUMERICAL RESULTS: SOLITONS 



A. Symmetric attraction-attraction system 
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FIG. 4: (Color online) Families of 2D solitons found in the symmetric attraction-attraction model 
(Ai = A2 = 1, Ai = A2 = 0). Symmetric, antisymmetric and asymmetric solutions are labeled by "s", 
"an" , and "as" , respectively. Dashed-dotted lines represent solutions (of all the three types) produced by 
the variational approximation, while solid and dotted lines depict, respectively, numerically found stable 
and unstable solutions, (a) The soliton's total power versus the chemical potential, (b) Ratio [see Eq. 
([9])] for the family of asymmetric solitons. The same conventions for labelling solution families of different 
types (variaitonal/numerical, stable/unstable, symmetric/antisymmetric/asymmetric) are adopted in other 

figures. 




FIG. 5: (Color online) Comparison of profiles of the 2D soliton in the attraction-attraction model, as found 

from the numerical solution and predicted by the variational approximation (solid and dashed curves, 
respectively). Shown are cross sections of the soliton along x- and 7/-axes, for \i = — 12. Norms of the two 

components of the soliton are N u « 4.5 and N v « 0.4. 



essentially different from what was observed in the ID system). The branch of the asymmetric solitons is 
also subject to the collapse, but this happens on a remote portion of the N(fjb) curve, which could not be 
shown in Fig. 0J 



B. Symmetric repulsion-repulsion system 



A generic example of families of solitons found in the RR model is shown in Fig. In subgap lb, the VA 
predicts the solutions very accurately. For this case, the comparison of numerical and variational asymmetric 
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(c) (d) 

FIG. 6: (Color online) Evolution of an antisymmetric 2D soliton (with \i = —16) in the 
attraction- attraction model, which is destabilized by the collapse, (a) The initial state at t = 0; (b) the 
formation of an asymmetric state at t = 4; (c) the collapse of the high- amplitude component and decay of 
the low-amplitude one, at t = 6; (d) the evolution of the entire pattern in the cross section along x = 0. 
Here and in similar figures presented below, panels display the side-by-side juxtaposition of the two 

components, \u\ and \v\. 

soliton profiles is displayed in Fig. [51 On the other hand, the VA also predicts asymmetric solutions in subgap 
2a, where such solutions could not be found in the numerical form, therefore the extension of the asymmetric 
branch into the latter subgap is an artifact of the variational method. It is noteworthy too that the stable 
branch of symmetric solitons crosses the Bloch band separating subgaps la and lb, thus providing for the 
first example of 2D embedded solitons in a self-defocusing medium. These embedded solitons are stable, as 
illustrated by an example of the evolution of a perturbed soliton shown in Fig. [9j 

The situation in the RR system is closer to what was found for its ID counterpart in Ref. [12] (because 
collapse does not occur with repulsive nonlinearity ) : the antisymmetric branch undergoes a supercritical 
bifurcation, which gives rise to stable asymmetric solitons, while the symmetric branch does not bifurcate 
and remains always stable (the VK criterion does not apply to gap solitons in models with the repulsive 
nonlinearity). Note that, like in Fig. [H the bifurcation point is very accurately predicted by the VA. The 
asymmetric solitons are stable whenever they exist, hence the bist ability between symmetric and antisymmet- 
ric or asymmetric solitons takes place here. The antisymmetric solitons, destabilized by the SSB bifurcation, 
transform themselves into persistent localized breathers, see an example in Fig. [10J The breather features 
oscillations in its two components, with equal amplitudes and phase shift tt/2 between them. 

Lastly, stable asymmetric solitons have also been found in higher bandgaps, an example of which is given 
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(a) (b) 

FIG. 7: (Color online) The same as in Fig. [H but for families of solitons in the symmetric 
repulsion-repulsion model (Ai = A2 = — 1, A = 0). The portion of the asymmetric-soliton branch in subgap 

2a is an artifact of the variational approximation. 
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FIG. 8: (Color online) Comparison of profiles of asymmetric solitons, as predicted by the VA and found in 
the numerical form (dashed and solid lines, respectively), in the symmetric repulsive-repulsive model. In 

this case, \i = —4, and N u w 10.5, N v « 0.6. 

in Fig. [TU Typically, the solitons in higher bandgaps are less tightly bound, featuring well-pronounced 
sidelobes. 



IV. NUMERICAL RESULTS: VORTICES 
A. Symmetric attraction-attraction system 

In the model based on the single-component GPE in two dimensions with the OL and attractive nonlin- 
earity, stable solutions in the form of localized vortices with topological charge ( "spin" ) S = 1 were reported 
in Refs. 28} and [33], and their (also stable) counterparts with S > 2 were found in Ref. [35[. Before 
proceeding to the new problem of the SSB in two-component vortices, it is relevant to recapitulate basic 
results obtained for vortex solitons in the single-component model. Figure [12] shows a generic example of 
families of vortex-soliton solutions with S = 1 in this model. The most compact "crater-shaped" vortices, 
which feature a single peak with a hole in the center, are always unstable. However, the vortices built as 4- 
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FIG. 9: (Color online) Evolution of a perturbed symmetric embedded soliton in the symmetric 
repulsion-repulsion model is shown in the cross-section along x = 0. The chemical potential of the 
unperturbed soliton is \i = —7.54, placing it into the Bloch zone separating subgaps la and lb. The norm 

of the soliton is N « 8. The soliton is obviously stable. 




FIG. 10: (Color online) Spontaneous transformation of an unstable antisymmetric soliton (with \i = —3.5) 
into a persistent breather, in the symmetric repulsion-repulsion model (Ai = A2 = — 1, A = 0). 




FIG. 11: (Color online) An example of a stable asymmetric soliton found, in the symmetric 
repulsion-repulsion system, in subgap 2a, with \i = 1. 
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(c) (d) 

FIG. 12: (Color online) Localized vortices with topological charge 5 = 1 in the single-component model 
with attraction (A = 1). (a) Three solution families (unstable single-peak crater-shaped vortices, and stable 
4- and 8-peak vortex complexes) are represented by the dependence of their norm on the chemical 
potential, (b)-(d) Examples of the three species of localized vortices, with (/i = —4.5, TV w 20), 
(/i = -20, N « 76), and (// = -12, N « 78), respectively. 



and 8-peak complexes, with phase shifts, respectively, tt/2 or 7r/4 between the peaks (which corresponds to 
the net phase circulation of 27r, i.e., S = 1), constitute entirely stable families (Fig. [T2l includes examples of 
all the three species of localized vortices). These complexes are quite similar to examples of stable vortices 
reported in Refs. [28[ and [33( | . 

We have found the SSB of two-component vortex solitons in the symmetric AA system (Ai = A2 = +1, 
A = 0). Families of 4- and 8-peak vortices of the symmetric, antisymmetric and asymmetric types, found 
in the coupled system, are shown in Figs. [13] and [HJ along with examples of respective stable asymmetric 
vortices. 

The stability of these solutions was identified, as above, by direct simulations. Symmetric vortices get 
destabilized by the bifurcation and tend to spontaneously rearrange into asymmetric ones, which emerge as 
stable solutions, while antisymmetric vortices do not bifurcate. At large values of the norm, the localized 
vortices are subject to collapse. 
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(b) 



FIG. 13: (Color online) 4-peak vortices with topological charge S = 1, in the symmetric 
attraction-attraction system, with k = 2. (a) Families of the symmetric, antisymmetric, and asymmetric 
vortex solutions represented by N(/i) curves, (b) An example of a stable asymmetric vortex, for \i = — 18 

and N = 36.5. 





(b) 



FIG. 14: (Color online) The same as in Fig. [131 but for 8-peak vortices. In panel (b), an example of the 
stable asymmetric vortex is displayed for ji = —18, N = 73.5. 



B. Symmetric repulsion-repulsion system 

In the repulsive model with the OL, localized vortices may be treated as a species of 2D solitons of the gap 
type [3, [34)]. First, in Fig. [15] we present families of vortex solitons in the single-component model. Similar to 
what was shown above for the attractive model, the single-peak (crater-shaped) vortices are unstable, while 
multiple-peak vortical complexes are stable. However, in the repulsive model only the (unstable) single-peak 
and (stable) 8-peak structures carry topological charge S = 1, while the 4-peak entity is a complex bound 
state of several vortices. Indeed, the phase distribution in the latter state, displayed in Fig. [16l suggests 
that it may be interpreted as a structure built of of five vortices: one, with S = 1, is located in the center, 
being surrounded by four constituent vortices, each carrying charge S = — 1. 

In the coupled (two-component) RR system, we observed SSB in both types of stable vortex states, 4- 
peak and 8-peak ones, as shown in Fig. [TT1 Similar to what was demonstrated above for the solitons, 
the asymmetric branch bifurcates from the antisymmetric one, while the family of the symmetric vortices 
does not bifurcate and remains entirely stable, thus giving rise to the bist ability, together with the stable 
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(b) 





(c) 



(d) 



FIG. 15: (Color online) Localized vortices in the single- component model with repulsion (A = —1). (a) 
Norm versus the chemical potential for three species of vortex solitons: single-peak, 4-peak, and 8-peak 
ones, (b)-(d): Examples of vortices of these types, with (/i = 1.5, N w 35), (fi = —2.5, N w 135), and 

(/I = —3.5, N w 208), respectively. 




(a) 



FIG. 16: (Color online) The phase pattern in the 4-peak vortex complex in the single-component repulsive 

model, whose density profile is displayed in Fig. [T5l c). 
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FIG. 18: (Color online) A stable asymmetric vortex with topological charge S = 2, found in the symmetric 

repulsion-repulsion system. In this case, \i = 1 (which falls in subgap 2a). Norms of the components are 
N u & 272, N v & 13. Note the difference in the x- and y- scales in the figure (in fact, the vortex is invariant 

with respect to x <-> y). 




asymmetric vortices. Antisymmetric two-component vortices destabilized by the SSB bifurcation develop 
into persistent breathers. A noteworthy feature of Fig. [171 is the fact that the branch of symmetric solutions 
(both 4- and 8-peak ones) crosses the Bloch band between subgaps la and lb, which provides for the first 
example of (stable) embedded vortex solitons. 

Stable asymmetric vortices residing in higher bandgaps, as well as stable vortices with higher vorticity, 
S > 1, have been found too. In particular, Fig. [18] displays an example of a stable 8-peaked vortex with 
S = 2, found in subgap 2a. 



V. ASYMMETRIC MODELS 



As explained in Introduction, the symmetry of coupled equations ([Taj) may be broken by opposite signs of 
the nonlinearity in the two cores (in the AR system, with Ai = — A2 = +1), or by mismatches between the 
lattices in them, Ai 5 2 7^ 0. We did not aim to perform an exhaustive analysis of the asymmetric models, but 
examples of stable solitons in these models have been found. 
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FIG. 19: (Color online) Stable solitons in the coupled attraction-repulsion system 
(Ai = — A2 = 1, Ai ? 2 = 0). The solitons' profiles are even functions of x and y, and are invariant with 
respect to transformation x <-> y. (a) A soliton with the dominant self-attractive component. It has 
ji w —13 (which falls in the semi-infinite gap), and N u w b.7,N v « 0.26. (b) A soliton with the dominant 
self-repulsive component, found with \i « —5.7 (which belongs to subgap la) and N u w 0.55, iV v « 5.35. 



Figure [19] shows stable solitons in the AR system. Solutions of two types have been found in it: one with 
a dominant self-repulsive component sitting in subgap la, and another one with the dominant self- attractive 
component found, as should be expected, in the semi-infinite gap. Stability of the solitons have been verified 
by direct simulations. 

Stable 8-peak vortices with S = 1 have also been found in the AR model. Vortices with dominant self- 
attractive and self-repulsive components are located in the semi-infinite gap and subgap la, respectively. 
Their examples are not shown here, as they are very similar to the 8-peak asymmetric solitons found in the 
A A and RR models, that were reported above. 

We have also studied the system with the mismatch between the OLs, concentrating on two most interesting 
cases, with Ai = 7r, A2 = or Ai = A2 = tt in Eqs. ([Taj) . In either case, the phase mismatch was given the 
largest possible value (tt). As in the ID model with nonzero mismatch [12|, quasi-symmetric and asymmetric 
solutions have been identified, as states with N u = N v and N u 7^ N v , respectively. In the mismatched 
system, asymmetric solitons have peaks in the two components located at the same position, while in quasi- 
symmetric solutions the peaks are slightly separated. Typical profiles of quasi- symmetric and asymmetric 
solutions in the misaligned AA system are shown in Fig. [20] for Ai = A2 = tt (this choice corresponds to the 
largest diagonal mismatch). Direct simulations demonstrate that the quasi-symmetric solitons are always 
unstable in this case (an explanation of this feature is given below), while the asymmetric solitons are stable, 
beneath the collapse threshold. 

Further, Fig. [2T1 depicts the dependence of the asymmetry ratio, defined as per Eq. ([9]), on coupling coef- 
ficient ft. In the case of the maximum horizontal mismatch, Ai = 7r, A2 = 0, we observe that the bifurcation 
which breaks the (quasi-) symmetry of solitons at a critical value of ft is replaced by the pseudo-bifurcation, 
i.e., the branch of asymmetric solutions asymptotically approaches its quasi-symmetric counterpart with 
the increase of ft, but the bifurcation does not happen, since the two branches never merge. A similar 
phenomenon was reported in the ID model with A = tt [l2[. The replacement of the bifurcation by the 
pseudo-bifurcation is observed only in the case of the largest mismatch, Ai = tt and/or A2 = n, and it ex- 
plains the above-mentioned total instability of the family of quasi-symmetric solitons, as the corresponding 
branch never has a chance to get stabilized by the inverse bifurcation (if one considers the evolution of the 
solutions with the increase of ft at fixed total norm N). 

A new feature specific to the 2D system is that, in the system with the largest diagonal mismatch, 
Ai = A2 = 7r, the numerical procedure suddenly ceases to converge at some critical value of ft, and no 
asymmetric solitons are found past this point. For instance, at ft = 15.485 a stable asymmetric soliton with 
a regular profile can be found, but when the coupling increases to ft = 15.490, the soliton is lost. The shape 
of the 9 (ft) dependence in Fig. 121T b) strongly suggests that the branch of asymmetric solitons terminates, 
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FIG. 20: (Color online) Soliton profiles in the attraction-attraction model with the maximum diagonal 
mismatch between the lattices in the cores, Ai = A2 = 7r, for n = 10, and TV = 5. (a) A stable asymmetric 

soliton; (b) an unstable quasi-symmetric soliton. 
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(a) (b) 

FIG. 21: (Color online) The asymmetry ratio, defined as per Eq. ([9]), as a function of the linear-coupling 
constant in the misaligned (Ai = 7r, A2 = 0) attraction-attraction model with fixed total norm of the 
two-component soliton, N = 5. (a) The pseudo-bifurcation, (b) The termination of the family of 

asymmetric soliton solutions. 

in this case, due to a tangent bifurcation, which results from collision and mutual annihilation of the branch 
with an additional branch of unstable solutions, that we did not aim to find. 

Other types of stable solitons and vortices have also been found in the misaligned system of the AA type. 
In particular, an example of a soliton with two peaks in one component and a single peak in the other is 
shown in Fig. [22j 

Stable asymmetric solitons in the RR system with misaligned lattices have been found too. A typical 
stable asymmetric soliton in this system is shown in Fig. [23] 

VI. CONCLUSION 



We have introduced experimentally relevant models of two stacked flat BEC traps, with attractive or 
repulsive interactions between atoms. Each flat trap carries an OL (optical lattice), stable 2D solitons being 
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FIG. 23: (Color online) An example of a strongly asymmetric stable soliton in the diagonally mismatched 
(Ai = A2 = 7r) repulsion-repulsion system with n = 10. This soliton has n = — 11, and 7V n = 4.4, 7V V = 1.6. 



impossible without it. First, it was demonstrated that the linear coupling splits every finite bandgap in 
the linear spectrum of the single-component model into two subgaps, which are separated by narrow Bloch 
bands. The main issue, addressed in this work by means of the VA (variational approximation) and numerical 
methods, is SSB (spontaneous symmetry breaking) in families of solitons and localized vortices, as a result 
of the competition of three factors: the attractive or repulsive nonlinearity in each core, the action of the 
OL potential, and the linear coupling between the cores. Similar to what was found in the zero-dimensional 
setting flQl l (t his means a double- well potential without any transverse dimension) and one-dimensional 
models [12|, [l9j, the SSB occurs in families of symmetric or antisymmetric states, in the case of the attractive 
or repulsive nonlinearity, respectively. In either case, the corresponding bifurcation destabilizes the branch 
of symmetric or antisymmetric solitons or vortices, giving rise to a stable branch of asymmetric solutions. 
The VA, based on the Gaussian ansatz, yields an accurate prediction for the bifurcation and the branch of 
asymmetric solitons generated by the bifurcation. A feature specific to the 2D setting is the termination of 
all stable branches in the AA (attraction-attraction) system due to the onset of collapse. 

Stable asymmetric solitons sitting in higher finite bandgaps, and localized vortices with multiple values of 
the topological charge have been found too. In addition, the models considered in this work give rise to first 
examples of (stable) embedded solitons and embedded vortices (localized states found inside Bloch bands 
separating the subgaps) in any 2D setting. 

Solitons and localized vortices were also considered in the linearly-coupled systems whose symmetry is 
broken by opposite signs of the nonlinearity in the two traps, or by a mismatch between the OLs in them. In 
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the former case, the system gives rise to two distinct types of stable solitons and vortices, which are dominated 
by either a self-attractive component or self-repulsive one, which sit, respectively, in the semi-infinite gap, 
or in a finite bandgap. In the system with mismatched lattices, the phenomenon of the pseudo-bifurcation 
(which was recently reported in the ID system [12[) was found: when the mismatch takes the largest value (tt) 
in any direction (horizontal or diagonal), the SSB bifurcation fails to happen, as branches of asymmetric and 
quasi- symmetric (or quasi- antisymmetric) solutions asymptotically approach each other, but never merge. 

This work was supported, in a part, by the Israel Science Foundation through the Center-of- Excellence 
Grant No. 8006/03 and the German-Israel Foundation (GIF) through Grant No. 149/2006. 
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